#include "2BWTSearch.h"
#include "2BWTSearch_cpu.h"
#include "mian_util.h"
#include "cuda_header.h"
#include <algorithm>
#include "coproc.h"
//#include "Index.h"
//#include "PausedQueryBuffer.h"


using namespace std;

char getComplement(char base) {

        switch(base) {
                case 'a':
                        return 't';
                case 'A':
                        return 'T';
                case 't':
                        return 'a';
                case 'T':
                        return 'A';
                case 'c':
                        return 'g';
                case 'C':
                        return 'G';
                case 'g':
                        return 'c';
                case 'G':
                        return 'C';
                default: {
			printf("character : %c \n", base);
                        ERROR_EXIT("unknown base");
		}
        }
}
void reverse_complement(char* outRead, const char* inRead, const int readLength) {
        for(int i = 0; i < readLength; i++) {
                outRead[i] = getComplement(inRead[i]);
        }
        reverse(outRead, outRead + readLength);
}



void copyFileName(char * name, char * path, char * suffix)
{
	strcat(name, path);
	strcat(name,suffix);
}

char decodeInt(int i)
{
	switch(i)
	{
		case(0):
			return 'A';
			break;
		case(1):
			return'C';
			break;
		case(2):
			return'G';
			break;
		case(3):
			return'T';
			break;
	}
}

char EncodeInt(char a)
{
	switch(a)
	{
	case'A':
		return 0;
	case'C':
		return 1;
	case'G':
		return 2;
	case'T':
		return 3;
	}
	return 0;
}


uint64 reportNumMatch(const long* numMatch, const int numQuery) {
	int numDiscard = 0;
	int numAlign = 0;
	int numQueryAligned = 0;
	
	for(int i = 0; i < numQuery; i++) {
		if(numMatch[i] == -1) {
			numDiscard++;
		} else {
			numAlign += numMatch[i];
			if(numMatch[i] > 0) {
				numQueryAligned ++;
			}
		}
	}
	
	printf("reported by reportNumMatch: numDiscard: %d\n", numDiscard);
	printf("reported by reportNumMatch: numQueryAligned: %d\n", numQueryAligned);
	printf("reported by reportNumMatch: numAlign: %d\n", numAlign);

	return numAlign;
}

void LoadIndex(BWT * bwt, BWT *revbwt, char* path)
{
	int i;
	FILE *bwtCodeFile,*saValueFile, *revbwtCodeFile, *occValueFile,*revoccValueFile, *gpuOccValueFile, *gpurevOccValueFile= NULL;
	unsigned int tmp;
	unsigned int bwtCodeLengthInFile;
	const int MaxFileNameLength = 100;
//	char path[] = "/home/lumian/data/human_ref/ref_good.fa.index";
//	char path[] = "/home/lumian/data/human_ref/2bwt-index/hg19.fa.index";
	//char path[] = "data/ref.fa.index";
	char bwtCodeFileNameSuffix[] = ".bwt";
	//	char saValueFileNameSuffix[] = "/ref.fa.index.sa";
	char revbwtCodeFileNameSuffix[] = ".rev.bwt";
	char lookupTableFileNameSuffix[] = ".lkt";
	char revLookupTableFileNameSuffix[] = ".rev.lkt";
	char occValueFileNameSuffix[] = ".fmv";
	char revoccValueFileNameSuffix[] = ".rev.fmv";
	char gpuoccValueFileNameSuffix[] = ".fmv.gpu";
	char gpurevoccValueFileNameSuffix[] = ".rev.fmv.gpu";
//	char gpuoccValueFileNameSuffix[] = ".fmv.gpu";
//	char gpurevoccValueFileNameSuffix[] = ".rev.fmv.gpu";


	char bwtCodeFileName[MaxFileNameLength]="";
	char saValueFileName[MaxFileNameLength]="";
	char revbwtCodeFileName[MaxFileNameLength]="";
	char lookupTableFileName[MaxFileNameLength]="";
	char revLookupTableFileName[MaxFileNameLength]="";
	char occValueFileName[MaxFileNameLength]="";
	char revoccValueFileName[MaxFileNameLength]="";

	//char gpuOccValueFileName[] = "/home/lumian/data/human_ref/hg19_OccValueFileGPU";

	char gpuOccValueFileName[MaxFileNameLength] = "";
	char gpurevOccValueFileName[MaxFileNameLength] = "";

        //char gpuOccValueFileName[] = "OccValueFileGPU";
        //char gpurevOccValueFileName[] = "RevOccValueFileGPU";

	copyFileName(bwtCodeFileName, path, bwtCodeFileNameSuffix);
	//	copyFileName(saValueFileName, path, saValueFileNameSuffix);
	copyFileName(revbwtCodeFileName, path, revbwtCodeFileNameSuffix);
	copyFileName(lookupTableFileName, path, lookupTableFileNameSuffix);
	copyFileName(revLookupTableFileName, path, revLookupTableFileNameSuffix);
	copyFileName(occValueFileName, path, occValueFileNameSuffix);
	copyFileName(revoccValueFileName, path, revoccValueFileNameSuffix);
	copyFileName(gpuOccValueFileName, path, gpuoccValueFileNameSuffix);
	copyFileName(gpurevOccValueFileName, path, gpurevoccValueFileNameSuffix);

	unsigned int *tempBwtCode;
	unsigned int * tempPacValue;

	bwtCodeFile = (FILE*)fopen(bwtCodeFileName, "rb");
	if (bwtCodeFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open bwtCodeFile%s!\n", bwtCodeFileName);
		exit(1);
	}

	revbwtCodeFile = fopen(revbwtCodeFileName, "rb");
	if (revbwtCodeFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open revbwtCodeFile!\n");
		exit(1);
	}
#if LAYER_INDEX_2
	occValueFile = (FILE*)fopen(occValueFileName, "rb");
	if (occValueFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open occFile!\n");
		exit(1);
	}

	revoccValueFile = fopen(revoccValueFileName, "rb");
	if (revoccValueFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open revoccFile!\n");
		exit(1);
	}
#endif

#if LAYER_INDEX_1
	gpuOccValueFile = (FILE *)fopen(gpuOccValueFileName, "rb");
	gpurevOccValueFile = (FILE *)fopen(gpurevOccValueFileName,"rb");
#endif
	/*
	   if (saValueFileName != NULL && saValueFileName[0] != '\0' && saValueFileName[0] != '-') {
	   saValueFile = (FILE*)fopen(saValueFileName, "rb");
	   if (saValueFile == NULL) {
	   fprintf(stderr, "BWTLoad() : cannot open saValueFile!\n");
	   fprintf(stderr, "saValueFileName:%s\n", saValueFileName);
	   exit(1);
	   }
	   }

*/
	printf("reference file: %s\n",path);
	bwt->countValue = (unsigned int *)malloc((ALPHABET_SIZE + 1)*sizeof(int));
	fread(&bwt->inverseSa0, sizeof(unsigned int), 1, bwtCodeFile);
	fread(&bwt->countValue[1], sizeof(unsigned int), ALPHABET_SIZE, bwtCodeFile);
	Initialize(bwt, bwt->countValue[ALPHABET_SIZE]);
	bwt->countValue[0] = 0;

	bwt->bwtSizeInWord = BWTResidentSizeInWord(bwt->textLength) + WORD_BETWEEN_OCC / 2;
	bwtCodeLengthInFile = BWTFileSizeInWord(bwt->textLength);
	tempBwtCode = (unsigned int *)malloc(bwt->bwtSizeInWord * sizeof(unsigned int));
	printf("reference length: %u\n", bwt->textLength);	
	

	fread(tempBwtCode, sizeof(unsigned int), bwtCodeLengthInFile, bwtCodeFile);
	fclose(bwtCodeFile);

	DecodeBWT(bwt, tempBwtCode);
	BWTClearTrailingBwtCode(bwt);

	LoadLookupTable(bwt->lookupTable, lookupTableFileName, TABLE_SIZE);

#if LAYER_INDEX_2
	for(i=0; i<ALPHABET_SIZE+1; i++)
		fread(&tmp,sizeof(unsigned int), 1, occValueFile);
	fread(bwt->occValue, sizeof(unsigned int), bwt->occSizeInWord, occValueFile);
	fread(bwt->occMajorValue, sizeof(unsigned int), bwt->occMajorSizeInWord, occValueFile);
	fclose(occValueFile);
#endif

#if LAYER_INDEX_1
	uint32 temp;
	fread(&temp, sizeof(unsigned int), 1, gpuOccValueFile);
	printf("sa0? : %u\n", temp);
	for(i=0; i<ALPHABET_SIZE; i++){
		fread(&temp, sizeof(unsigned int), 1, gpuOccValueFile);
		printf("%u\n", temp);
	}
	unsigned int gpuOccValueSize = (bwt->textLength + GPU_OCC_INTERVAL - 1) / GPU_OCC_INTERVAL * ALPHABET_SIZE;	
	bwt->GPUOccValue = (unsigned int *)malloc(gpuOccValueSize * sizeof(int));
	fread(bwt->GPUOccValue, sizeof(int), gpuOccValueSize, gpuOccValueFile);
	fclose(gpuOccValueFile);
	printf("interval? : %u\n", bwt->GPUOccValue[0]);
	bwt->GPUOccValueInWord = gpuOccValueSize;

#endif
	
	revbwt->countValue = (unsigned int *)malloc((ALPHABET_SIZE + 1)*sizeof(int));
	fread(&revbwt->inverseSa0, sizeof(unsigned int), 1, revbwtCodeFile);
	fread(&revbwt->countValue[1], sizeof(unsigned int), ALPHABET_SIZE, revbwtCodeFile);
	Initialize(revbwt, revbwt->countValue[ALPHABET_SIZE]);
	revbwt->countValue[0] = 0;

	bwtCodeLengthInFile = BWTFileSizeInWord(revbwt->textLength);
	revbwt->bwtSizeInWord = BWTResidentSizeInWord(revbwt->textLength) + WORD_BETWEEN_OCC / 2;
	tempBwtCode = (unsigned int *)malloc(bwt->bwtSizeInWord * sizeof(unsigned int));
	fread(tempBwtCode, sizeof(unsigned int), bwtCodeLengthInFile, revbwtCodeFile);
	fclose(revbwtCodeFile);

	DecodeBWT(revbwt, tempBwtCode);
	BWTClearTrailingBwtCode(revbwt);
	LoadLookupTable(revbwt->lookupTable,revLookupTableFileName, TABLE_SIZE);

#if  LAYER_INDEX_1	
	fread(&temp, sizeof(unsigned int), 1, gpurevOccValueFile);
	for(i=0; i<ALPHABET_SIZE; i++)
		fread(&temp, sizeof(unsigned int), 1, gpurevOccValueFile);
	revbwt->GPUOccValue = (unsigned int *)malloc(gpuOccValueSize * sizeof(int));
	fread(revbwt->GPUOccValue, sizeof(int), gpuOccValueSize, gpurevOccValueFile);
	fclose(gpurevOccValueFile);
	revbwt->GPUOccValueInWord = gpuOccValueSize;
#endif
	
#if LAYER_INDEX_2
	for(i=0; i<ALPHABET_SIZE+1; i++)
		fread(&tmp,sizeof(unsigned int), 1,revoccValueFile);
	fread(revbwt->occValue, sizeof(unsigned int), revbwt->occSizeInWord, revoccValueFile);
	fread(revbwt->occMajorValue, sizeof(unsigned int), revbwt->occMajorSizeInWord, revoccValueFile);
	fclose(revoccValueFile);
#endif
}

void parseInputFile(char * inputFileName, int ReadLength, Query *Reads, 
					unsigned int* ReadNum, unsigned int * Queries, 
					const bool reverse = false)
{
	char c;
	unsigned int i=0, j=0, k=0, pos = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = QUERY_LENGTH + 10;
	char ReadName[maxReadNameLength];
	unsigned int querySizeInWord = getQuerySizeInWord(QUERY_LENGTH);
	FILE * fin = fopen(inputFileName, "r");
	char temp[200];
	bool FASTQ = false;
	printf("queryFile: %s\n", inputFileName);
	if(strstr(inputFileName, "fastq") != NULL){
		FASTQ = true;
		printf("fastq\n");
	}
	char* revQuery = (char*)malloc(sizeof(char)*ReadLength);
	
	while(i<QUERY_NUM && !feof(fin))
	{
		c = fgetc(fin);
		if(FASTQ){
		while(c!='@' && !feof(fin)){
			c = fgetc(fin);
		}}
		else{
		while(c!='>' && !feof(fin)){
			c = fgetc(fin);
		}}
		if(feof(fin))
			break;
		Reads[i].queryName = new char[maxReadNameLength];
		fgets(Reads[i].queryName,maxReadNameLength,fin);
		for(j=0; j<maxReadNameLength; j++){
			if(Reads[i].queryName[j] == '\n'){
				Reads[i].queryName[j] = '\0';
				break;
			}
		}
		Reads[i].query = new char[maxReadLength];
		pos = 0;
		fgets((char *)Reads[i].query, maxReadLength, fin);
		for(j=0; j<ReadLength; j++){
			c = Reads[i].query[j];
			if(c == 'a' || c =='c' || c=='g' || c=='t'){
				Reads[i].query[pos++] = toupper(c);
			}
			else if( c == 'A' || c =='C' || c=='G' || c=='T')
				Reads[i].query[pos++] = c;
			else if(c == 'N' || c=='n')
				Reads[i].query[pos++] = 'G';

		}
		Reads[i].query[pos] = '\0';	
		
		
		if(reverse == true) {
			reverse_complement(revQuery, Reads[i].query, ReadLength);
			memcpy(Reads[i].query, revQuery, sizeof(char)*ReadLength);
		}
		
		
//		printf("%s\n",Reads[i].query);	
		//pack Reads[i].query to Queries 
		pos = 0;
                for(j=0; j<querySizeInWord; j++) {
                        for(k=0; k<CHAR_PER_WORD && pos<QUERY_LENGTH; k++) {
                                Queries[i * querySizeInWord + j] <<= 2;
                                Queries[i * querySizeInWord + j] |= (EncodeInt(Reads[i].query[pos]) & 0x3);
                                pos++;
                        }
                        if(pos == QUERY_LENGTH) {
                                for(; k<CHAR_PER_WORD; k++)
                                        Queries[i * querySizeInWord + j] <<= 2;
                        }
                }
		if(FASTQ) {
			fgets(temp,200, fin);
			fgets(temp,200, fin);
		}
		i++;
	}
	*ReadNum = i;
	fclose(fin);
	
	free(revQuery);
}

uint64 getNumAlign(uint3* sa, const uint64 numSA) {
	uint64 numAlign = 0;
	for(uint64 i = 0; i < numSA; i++) {
		numAlign += (sa[i].z - sa[i].y + 1);
	}
	
	return numAlign;
}

void sortQuery(uint3* cpuSA, unsigned int* Queries, const int numQuery, 
				BWT* bwt, BWT* revbwt) {
	preprocess(cpuSA, Queries, numQuery, bwt, revbwt);
	
	FILE* file = NULL;
	safe_fopen(&file, "query.bin", "wb");
	safe_fwrite(Queries, sizeof(unsigned int), QUERY_SIZE_IN_WORD*numQuery, file);
	safe_fclose(file);
}


void parseBinInputFile(const char* fileName, unsigned int* Queries, const int numQuery) {
	FILE* file = NULL;
	safe_fopen(&file, fileName, "rb");
	safe_fread(Queries, sizeof(unsigned int), QUERY_SIZE_IN_WORD*numQuery, file);
	safe_fclose(file);
}


void calculateNumMatch(long* numMatch, const uint3* sa, const uint64 numSA) {
	for(int i = 0; i < numSA; i++) {
		numMatch[sa[i].x] += (sa[i].z - sa[i].y + 1);
	}
}


int main(int argc, char* argv[])
{

#if LAYER_INDEX_1      
	printf("LAYER_INDEX_1\n"); 
#endif

#if LAYER_INDEX_2
	printf("LAYER_INDEX_2\n");
#endif

	int i,j;
	int matchCounter = 0;
	unsigned int matchRead = 0;
	Query* Reads = new Query[QUERY_NUM];
	unsigned int ReadNum = 0;
	unsigned int timer = 0;
	char* indexFileName = argv[1];
	char* queryFileName = argv[2];

	BWT *bwt, *revbwt;
	bwt = (BWT *)malloc(sizeof(BWT));
	bwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	revbwt = (BWT *)malloc(sizeof(BWT));
	revbwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	LoadIndex(bwt, revbwt, indexFileName);
	printf("Load index done. \n");
	
	unsigned int querySizeInWord = getQuerySizeInWord(QUERY_LENGTH);
	//both positive and negative strands, thus 2x
	unsigned int * Queries = 
			(unsigned int *)malloc(2*QUERY_NUM * querySizeInWord * sizeof(int));

	uint64 numSA1, numSA2;
	uint64 numAlign1, numAlign2;


	//char* queryFileName = "/home/luo_group/data/human_reads/SRR001113_1M.fasta";	//35 bp
	//char* queryFileName = "/home/lumian/data/human_reads/SRR003092_1M.fastq";	//50 bp
	//char* queryFileName = "/home/lumian/data/human_reads/SRR073318_1M.fasta";	//67 bp
	//char* queryFileName = "/home/lumian/data/human_reads/SRR003084_1M_new.fasta";	//36 bp
	//char* queryFileName = "/home/lumian/data/human_reads/korea_1M.fasta"; //100 bp
	//read both positive and negative strand queries from the raw data
	printf("query file: %s\n", queryFileName);
	printf("input #queries: %d\n", QUERY_NUM);
	printf("query length: %d\n", QUERY_LENGTH);
	printf("query size in word: %d\n", QUERY_SIZE_IN_WORD);
	printf("MAX_NUM_ALIGN_PER_QUERY: %d\n", MAX_NUM_ALIGN_PER_QUERY);
	
	parseInputFile(queryFileName, QUERY_LENGTH, Reads, &ReadNum, 
					Queries, false);
	parseInputFile(queryFileName, QUERY_LENGTH, Reads, &ReadNum, 
					Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, true);
	/////////////////////////////////////////////////////////////////////////////
	//GPU-CPU coprocessing				
//	const int numCPUThread = 16;	
	//=============> find all valid matches
	//coproc2BWTSearch(Queries, QUERY_NUM*2, bwt, revbwt, numCPUThread);	
//	coproc2BWTSearch_simple(Queries, QUERY_NUM*2, bwt, revbwt, numCPUThread);	
	//coproc2BWTSearch_simple(Queries, QUERY_NUM, bwt, revbwt, numCPUThread);
	//coproc2BWTSearch_simple(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, bwt, revbwt, numCPUThread);

	//===================> with an threshold
/*	coproc2BWTSearchLimit_simple(Queries, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM,
                                     bwt, revbwt, numCPUThread);	
*/	
	////////////////////////////////////////////////////////////////////////////
	//without CPU-GPU coprocessing

	////////////////////////
	//initialization
	//CPU thread initialization
/*
	CPU2BWTInit(bwt, revbwt);
	uint64 cpuBufSize = 4*((uint64)(1024*1024*1024));
	uint3* cpuSA = (uint3*)malloc(cpuBufSize);
*/	
	//GPU thread initialization
	BWT * d_bwt = NULL;
	BWT * d_revbwt = NULL;	
	st_MemoryBuf* GPUMemBuf = NULL;
	const uint64 gpuBufSize = (2.5)*((uint64)(1024*1024*1024));
	GPU2BWTInit(&d_bwt, &d_revbwt, bwt, revbwt, &GPUMemBuf, gpuBufSize);

				

	//processing either by GPU or CPU
	//GPU===>
	timer = 0;
	startTimer(&timer);
//	#ifdef GPU_NO_THRESHOLD 
	//======> find all valid alignments
	uint64 nresult = 0;
	nresult += GPU2BWTSearch_shareBuf(Queries, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//	nresult += GPU2BWTSearch(Queries, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//	nresult += GPU2BWTSearch_noConcurent(Queries, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//	nresult += GPU2BWTSearch1Kernel(Queries, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
        GPUMemBuf->reset();
   	nresult += GPU2BWTSearch_shareBuf(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//  	nresult += GPU2BWTSearch(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//   	nresult += GPU2BWTSearch_noConcurent(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//   	nresult += GPU2BWTSearch1Kernel(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
        GPUMemBuf->reset();
	cudaThreadSynchronize();
	printf("\n\n===============> #alignment: %lu\n", nresult);
	endTimer(&timer, "GPU processing");


	//=======> with a threshold
/*	long* d_numAlign = NULL;
	GPUMALLOC((void**)&d_numAlign, sizeof(long)*QUERY_NUM);	
	cutilSafeCall(cudaMemset(d_numAlign, 0, sizeof(long)*QUERY_NUM));
        GPU2BWTSearchTwoStrandLimit(d_numAlign, Queries, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
        cudaThreadSynchronize();
	GPUMemBuf->close();
	long* gpu_numAlign = (long*)malloc(sizeof(long)*QUERY_NUM);
	FROMGPU(gpu_numAlign, d_numAlign, sizeof(long)*QUERY_NUM);
	reportNumMatch(gpu_numAlign, QUERY_NUM);
*/

	//CPU===>
/*	long* numMatch = (long*)malloc(sizeof(long)*QUERY_NUM);
	memset(numMatch, 0, sizeof(long)*QUERY_NUM);
	//======>find all valid
	timer = 0;
	startTimer(&timer);

	numSA1 = CPU2BWTSearch(cpuSA, Queries, QUERY_NUM, bwt, revbwt);
	numAlign1 = reportCPUResult(cpuSA, numSA1);
	numSA2 = CPU2BWTSearch(cpuSA, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, bwt, revbwt);
	numAlign2 = reportCPUResult(cpuSA, numSA2);
	printf("\n\n===============> #alignment (all valid hits): %lu\n", numAlign1 + numAlign2);
	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer, "CPU all matches");

	printf("\n\n\n");
*/
/*	//======>with the #match threshold
	printTitle("CPU2BWTSearchTwoStrandLimit");
	timer = 0;
	startTimer(&timer);
	CPU2BWTSearchTwoStrandLimit(numMatch, cpuSA, Queries, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM,
                                        QUERY_NUM, bwt, revbwt);
	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer, "CPU 2BWTSearchLimit");
	reportNumMatch(numMatch, QUERY_NUM);
*/

/*
	//===>check result for the limit implmentation
	long* numMatch_gold = (long*)malloc(sizeof(long)*QUERY_NUM);
	memset(numMatch_gold, 0, sizeof(long)*QUERY_NUM);
        numSA1 = CPU2BWTSearch(cpuSA, Queries, QUERY_NUM, bwt, revbwt);
	calculateNumMatch(numMatch_gold, cpuSA, numSA1);	
        numSA2 = CPU2BWTSearch(cpuSA, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, bwt, revbwt);
	calculateNumMatch(numMatch_gold, cpuSA, numSA2);
	uint64 gold_numAlign = 0;
	for(int i = 0; i < QUERY_NUM; i++) {
		if(numMatch_gold[i] <= MAX_NUM_ALIGN_PER_QUERY) {
			gold_numAlign += numMatch_gold[i];
		}
	}
	printf("numMatch_gold with the threshold %d: %lu\n", MAX_NUM_ALIGN_PER_QUERY, gold_numAlign);

*/

	//////////////////////////////////////////////////////////////////////////////
	

	//GPUBWT->close();
	free(bwt->lookupTable->table);
	free(revbwt->lookupTable->table);
	free(bwt);
	free(revbwt);
	free(Queries);


	return 0;
}

